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Adaptive  Control  of  Bivalirudin  in  the  Cardiac 

Intensive  Care  Unit* * * § 

Qi  Zhao,'  Thomas  Edrich,*  Member,  IEEE,  and  Ioannis  Ch.  Paschal  idis,-<  Fellow,  IEEE, 


Abstract — Bivalirudin  is  a  direct  thrombin  inhibitor  used  in  the 
cardiac  intensive  care  unit  when  heparin  is  contraindicated  due 
to  heparin-induced  thrombocytopenia.  Since  it  is  not  a  commonly 
used  drug,  clinical  experience  with  its  dosing  is  sparse.  In  earlier 
work  ([1])  we  developed  a  dynamic  system  model  that  accurately 
predicts  the  effect  of  bivalirudin  given  dosage  over  time  and  pa¬ 
tient  physiological  characteristics.  This  paper  develops  adaptive 
dosage  controllers  that  regulate  its  effect  to  desired  levels.  To 
that  end,  and  in  the  case  that  bivalirudin  model  parameters  are 
available,  we  develop  a  Model  Reference  Control  law.  In  the 
case  that  model  parameters  are  unknown,  an  indirect  Model 
Reference  Adaptive  Control  scheme  is  applied  to  estimate  model 
parameters  first  and  then  adapt  the  controller.  Alternatively, 
direct  Model  Reference  Adaptive  Control  is  applied  to  adapt 
the  controller  directly  without  estimating  model  parameters  first. 
Our  algorithms  are  validated  using  actual  patient  data  from  a 
large  hospital  in  the  Boston  area. 

Index  Terms — Bivalirudin,  Pharmacokinetics,  Parameter  Iden¬ 
tification,  Adaptive  Control. 


I.  Introduction 

The  US  health  care  system  is  viewed  as  costly  and  highly 
inefficient.  Among  the  many  reform  efforts,  the  meaningful 
use  of  Electronic  Health  Records  (EHRs)  is  invariably  seen  as 
a  key  to  improving  efficiency.  In  the  hospital,  the  digitization 
of  data  from  medical  devices  enables  the  development  of 
algorithms  that  can  automate  decision  making  and  facilitate 
treatment.  This  is  exactly  the  goal  of  this  paper  which  fo¬ 
cuses  on  automating  dosage  decisions  for  a  particular  drug 
-bivalirudin-  used  in  the  cardiac  Intensive  Care  Unit  (ICU). 

Bivalirudin  antagonizes  the  effect  of  thrombin  in  the  blood 
clotting  cascade,  thereby  preventing  complications  from  blood 
clotting.  It  is  currently  FDA-approved  for  short-term  anti¬ 
coagulation  of  patients  undergoing  cardiac  catheterization  to 
prevent  complications  due  to  undesired  blood  clots  [2],  [3], 
[4],  [5].  Bivalirudin  is  administered  to  patients  who  have  a 
contraindication  to  heparin.  It  is  infused  continuously,  and  is 
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GM093147,  by  the  NSF  under  grants  CNS-1239021  and  IIS-1237022,  by  the 
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eliminated  via  the  kidney  and  by  plasma  protease-metabolism. 
It  affects  the  coagulation  parameters  Partial  Thromboplastin 
Time  ( PTT )  and  the  International  Normalized  Ratio  (INR) 
in  a  dose-dependent  fashion.  The  PTT  value  is  measured  in 
seconds  and  it  will  be  used  as  the  output  one  wishes  to  regulate 
within  a  specific  range. 

Although  not  commonly  used  overall,  bivalirudin  is  finding 
increasing  use  in  the  ICU.  Residents  adjusting  the  infusion  rate 
of  bivalirudin  may  have  limited  experience,  thus,  risking  over- 
or  under-dosing.  Currently,  the  drug  is  regulated  empirically  or 
with  a  very  simple  nomogram  [6].  Adequate  anticoagulation  is 
necessary  to  avoid  the  risk  of  clot  formation,  but  overshooting 
increases  the  risk  of  bleeding.  Complicating  matters,  there 
is  considerable  inter-  and  intra-individual  variability  in  the 
response  to  bivalirudin.  Motivated  by  these  challenges,  in 
earlier  work  [7],  [8],  [1],  we  developed  methods  for  predicting 
future  PTT  values  given  past  infusion  rates  and  the  patient’s 
renal  and  liver  function  characteristics.  Related  work  has  used 
pharmacokinetic-pharmacodynamic  models  to  model  the  effect 
of  various  drugs,  see,  e.g.,  [9],  [10].  One  of  our  methods  in  [1] 
proposes  an  explicit  dynamic  system  model  which  was  shown 
to  produce  quite  accurate  results  when  tested  against  actual 
patient  data. 

In  this  paper,  we  pursue  what  we  view  as  the  natural  next 
step.  Leveraging  the  dynamic  system  model  from  [8],  [1],  we 
seek  to  synthesize  controllers  that  can  regulate  the  infusion  rate 
to  drive  PTT  within  a  desirable  range.  Other  methodologies 
such  as  expert  systems  have  also  been  used  for  controlling 
some  drugs  [11].  We  develop  two  types  of  control  laws.  First, 
assuming  that  a  dynamic  system  model  that  can  predict  PTT 
given  dosage  is  completely  characterized,  we  develop  a  Model 
Reference  Control  (MRC)  law.  Model  parameters,  however, 
may  be  viewed  as  not  known  with  certainty,  which  is  due 
to  modeling  errors  and  inter-  or  intra-individual  variability. 
To  overcome  this  problem,  we  develop  an  indirect  Model 
Reference  Adaptive  Control  (MRAC)  law  that  identifies  the 
model  parameters  first  and  then  adapts  the  controller  in  real¬ 
time.  Furthermore,  we  develop  a  direct  Model  Reference 
Adaptive  Control  law  that  adapts  the  controller  directly  with¬ 
out  estimating  model  parameters  first,  which  is  more  efficient. 
For  each  case,  we  present  analytical  and  numerical  evidence 
showing  that  the  controllers  do  drive  PTT  to  the  desirable 
range.  Our  numerical  validation  is  in  fact  done  using  actual 
patient  data  from  the  Brigham  and  Women’s  Hospital  -  a  large 
hospital  in  the  Boston  area. 

The  remainder  of  the  paper  is  organized  as  follows.  Sec.  II 
presents  the  dynamic  system  model  that  predicts  the  effect  of 
bivalirudin  given  dosage  and  patient  physiological  information. 
Sec.  Ill  presents  the  proposed  control  schemes;  Sec.  III-A  de- 
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velops  the  MRC  law  whereas  Sec.  III-B  develops  the  indirect 
MRAC  law  based  on  the  patient  model  but  with  unknown 
parameters.  Sec.  III-C  develops  the  direct  MRAC,  which  is 
more  efficient  in  adapting  the  controller.  Finally,  concluding 
remarks  appear  in  Sec.  IV. 

Notation:  We  use  bold  letters  to  denote  vectors  and  ma¬ 
trices;  typically  vectors  are  denoted  by  lower  case  letters 
and  matrices  by  upper  case  letters.  Vectors  are  assumed  to 
be  column  vectors  unless  explicitly  stated  otherwise.  For 
economy  of  space  we  write  x  =  (xi,  ■  ■  •  ,xn)  for  the  column 
vector  x  £  R".  In  addition,  we  use  lower  case  letters  to  denote 
time  domain  functions  (e.g.,  /(f)),  and  upper  case  letters  to 
denote  Laplace  transforms  (e.g.,  F(s)). 

II.  Dynamic  System  Model  Formulation 
A.  The  Model 

This  section  presents  a  Multiple  Input  Single  Output  (MISO) 
dynamic  system  model  that  attempts  to  explicitly  account  for 
the  way  bivalirudin  affects  PTT  in  patients.  The  model  was 
developed  and  validated  in  [1];  it  is  presented  here  briefly 
to  establish  the  notation  and  to  set  the  stage  for  the  control 
schemes  of  Sec.  III. 


pi2:initial  condition 


Fig.  1.  The  Multiple  Inputs  Single  Output  (MISO)  dynamic  system  model. 

The  key  quantity  (response)  we  would  like  to  predict  is  the 
PTT  at  each  time  t.  The  dynamic  model  structure  is  shown 
in  Fig  1.  There  are  9  inputs  which  are  denoted  by  ufit),  i  = 
1, . . .  ,9  and  correspond  to  important  physiological  variables 
used  as  predictors.  More  specifically,  inputs  . . .  ,ug(t) 

respectively  correspond  to: 

1)  Bival  rate  (mg/kg/h):  the  weight-based  bivalirudin  in¬ 
jection  rate. 

2)  GFR  (mL/min):  the  glomerular  filtration  rate. 

3)  SGOT  (Units/L):  the  Semm  Glutamic  Oxaloacetic 
Transaminase. 

4)  SGPT  (Units/L):  the  Serum  Glutamic  Pyruvic  Transam¬ 
inase. 

5)  TBILI  (mg/dL):  total  bilirubin. 

6)  ALB  (g/L):  Albumin. 


7)  PLT  (K/mcL):  Platelet  count. 

8)  HCT  (%):  Hematocrit. 

9)  FIB  (mg/dL):  Fibrinogen. 

More  detailed  description  of  these  physiological  variables  can 
be  found  in  [1], 

The  model  of  Fig.  1  has  a  single  output  -the  PTT  value- 
which  is  denoted  by  y(t).  There  is  also  a  single  state  variable 
denoted  by  x(t).  Overall  there  are  12  unknown  parameters: 
11  of  which  correspond  to  the  various  gains  and  are  denoted 
by  [3i,  i  =  1, . . . ,  11.  The  initial  condition  of  the  system  is  the 
12th  unknown  parameter  and  is  denoted  by  :c(0)  (/?i2 ) ■  The 
system  dynamics  are: 

x(t)  =  Ax(t)  +  Bu(f),  (1) 

y(t)  =  Cx(t)  +  Du(f), 

where  A  =  — /?3,  B  =  [/3\  0  •  •  •  0],  C  =  /32,  and 
D  =  [0  /?4  •  •  •  /3n].  Clearly,  this  is  a  Linear  Time  Invariant 
(LTI)  dynamic  system.  The  challenge  is  that  we  do  not  know 
the  model  parameters  and  we  only  have  non-uniform  sampled 
inputs  u(t),  and  clinical  observation  values  y{t)  at  certain 
times  t  for  each  patient.  It  is  therefore  necessary  to  translate  the 
continuous-time  system  dynamics  to  discrete-time  dynamics 
before  proceeding  with  parameter  identification. 

B.  Parameter  Identification 

Given  the  highly  non-uniform  sampled  data,  two  methods 
were  introduced  to  identify  model  parameters  in  [1],  First, 
after  converting  to  discrete-time  dynamics,  we  formulated  the 
parameter  identification  problem  as  the  nonlinear  optimization 
problem  of  minimizing  some  metric  of  fitness  to  a  training  set 
of  sampled  data. 

The  data  set  we  used  comes  from  the  STAR  (Surgical 
ICU  Translational  Research)  Center  at  Brigham  and  Women’s 
Hospital  in  Boston.  It  consists  of  records  for  233  patients 
including  the  predictors  and  the  output  PTT  value  sampled 
(non-uniformly)  over  time.  We  randomly  split  our  data  set  into 
a  training  set  corresponding  to  2/3  of  the  total  (155  patients) 
and  a  test  set  corresponding  to  1/3  of  the  total  (78  patients). 
We  use  the  former  to  identify  the  unknown  system  parameters 
and  the  latter  to  evaluate  the  performance  of  the  various  control 
laws  we  will  develop  in  subsequent  sections. 

More  specifically,  let  us  use  a  subscript  j  to  denote  the 
model  primitives,  i.e.,  the  state  Xj(t),  output  yj(t),  and  inputs 
Uj(t)  for  each  patient  j  =  1 .....  A’,  where  N  denotes  the 
number  of  patients  in  the  training  set.  To  distinguish  between 
measurements  of  'yj(t)  and  predictions  based  on  the  system 
dynamics  we  use  yj  (f )  for  the  former  and  y:j  (/)  for  the  latter. 
Suppose  for  each  patient  j  we  have  7j  measurements  at  times 
tp . . .  ,/J3,  where  we  adopt  the  convention  t®  =  0  for  all  j. 
We  can  then  formulate  a  nonlinear  optimization  problem  of 
minimizing  the  least  squares  error 

N  V* 
i=!  t=t) 
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Fig.  2.  In  this  dynamic  model,  the  bivalirudin  infusion  rate  u(t)  is  the  only  controllable  input.  d(t)  is  the  linear  combination  of  the  rest  of  the  inputs. 


Reference  Model 


Fig.  3.  Model  Reference  Adaptive  Control  (MRAC)  structure. 


subject  to  the  discretized  version  of  the  system  dynamics.  The 
decision  variables  are  x(0)(=  x(tj)  for  all  j )  and  the  parame¬ 
ters  Pi,  i  =  1, . . . ,  11.  One  can  easily  substitute  the  expressions 
from  the  constraints  (the  system  dynamics)  into  the  objective 
function  and  obtain  a  nonlinear  optimization  problem  with  no 
constraints  other  than  some  bounds  on  the  decision  variables. 
This  problem,  however,  is  non-convex.  We  applied  a  Quasi- 
Newton  method  (BFGS)  [12]  to  obtain  optimal  values  for  the 
model  parameters.  This  yielded  a  population-wide  model  in  the 
sense  that  its  parameters  produced  the  best  fit  with  the  sampled 
data.  Furthermore,  and  to  accommodate  variability  across 
patients,  we  used  a  recursive  estimation  method  (Extended 
Kalman  Filter)  to  estimate  the  parameter  values  that  best  fit  a 
given  individual  patient  in  real-time. 

III.  Bivalirudin  Control  Scheme 

We  now  turn  to  our  primary  goal  of  devising  a  proper 
controller  to  keep  the  PTT  value  in  the  range  of  60s-80s. 
According  to  clinical  experience,  this  range  is  safe  and  optimal 
for  cardiac  surgery  patients.  For  the  system  we  have  defined  in 
Eq.  (1),  note  that  the  only  controllable  part  is  the  bivalirudin 
infusion  rate.  The  rest  of  the  inputs  are  indicators  of  patients’ 
liver  (SGOT,  SGPT,  TBILI,  ALB)  and  renal  function  (GFR) 


and  include  some  metrics  related  to  the  blood  (PLT,  HCT, 
FIB).  Arguably,  these  variables  are  not  immediately  affected 
by  the  drug  but  change  over  a  longer  time-scale  than  the 
one  we  focus  on  for  controlling  PTT  through  the  infusion  of 
bivalirudin.  Therefore,  we  consider  them  as  non-controllable 
and  aggregate  them  in  a  variable  d(t)  which  is  their  linear  com¬ 
bination  with  the  appropriate  gains  P±, . . . ,  pu  (see  Fig.  2). 

Ideally,  we  want  to  design  a  reference  model  which  can 
generate  sufficient  but  safe  PTT  values  driven  by  a  reference 
input  signal.  Based  on  the  output  of  the  reference  model, 
we  want  to  drive  our  system  to  perform  similarly  to  the 
reference  model  by  a  proper  control  signal.  Motivated  by 
this,  we  adopt  the  so  called  continuous-time  MRAC  scheme. 
Fig.  3  shows  the  general  structure  of  this  control  law.  IF,,,  (s) 
denotes  an  ideal  reference  transfer  function  that  can  generate 
the  desired  reference  output  signal.  The  controllable  system  is 
represented  by  Gp(s ,  0*),  where  0*  is  a  parameter  vector.  The 
objective  is  to  design  a  controller  C{s,0*e),  parameterized  by 
Q*e,  to  generate  the  proper  control  signals  that  can  drive  the 
controllable  system  to  track  the  reference  output  values. 

Our  first  controller  is  an  MRC  law  that  is  designed  assuming 
that  the  system  parameters  6*  are  known. 
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A.  Model  Reference  Control  (MRC) 

By  observing  the  system  in  Fig.  2,  we  can  rewrite  the 
dynamics  of  a  particular  patient  as 

xp(t)  =  ~P3xp(t)  +  P\u(t),  (2) 

yP(t)  =  p2Xp(t)  +  d(t),  (3) 


||x(f;f0,xo)  —  xe ||  <  ee  t°\  Vf  >  to  whenever  ||x0  — 
xe||  <  5(e). 

Theorem  III.l  If  we  choose  am  >  0 ,  bm  0,  and  r(t)  =  Cr 
(constant),  the  reference  model  equilibrium  state  yme  = 
is  exponentially  stable. 


where  we  use  u(t),xp(t),yp(t)  to  denote  the  input  signal 
(bivalirudin  infusion  rate),  the  state  variable,  and  the  output 
signal  (PTT),  respectively,  and  where  d(t)  =  Pi+2Ui(t). 
Clinically,  since  the  renal/liver  functions  and  blood  metrics 
of  patients  do  not  vary  much  within  a  certain  period,  we  do 
not  need  to  measure  these  physiological  variables  continuously 
and  we  assume  that  they  are  constant  within  the  sample 
interval.  By  observing  the  clinical  data,  we  find  that  d(f)  is 
a  step-wise  signal.  Therefore,  we  assume  that  the  first  order 
derivative  of  d(t)  ( d(t ))  is  0  within  the  sample  interval.  By 
taking  the  derivative  on  both  sides  of  (3),  using  (2)  to  substitute 
for  xp(t),  and  using  (3)  to  eliminate  xp(t ),  we  obtain: 


yP(t)  =  ~P3  yP(t)  +  P1P2  u(t)  +  (3  3d(t).  (4) 


In  the  frequency  domain,  we  have 


TO 


P1P2  U(s)  +  P3D(s) 
s  +  (3 3 


where  the  system  output  yp(t)  (Yp(s)),  the  input  u(t)  ( U(s )), 
and  d(t)  ( D(s ))  can  be  observed.  Hence,  in  our  setting,  the 
system  transfer  function  is  Gp(s,0*)  =  Yp(s)/U(s)  and  it  is 
parameterized  by  (3\,  P2  and  (33.  In  this  case,  9*p  =  ((3i,  P2, 
(3s)’ 

Next,  we  design  a  reference  transfer  function  Wm(s).  We 
take  Wm(s)  to  be  a  first-order  LTI  system  driven  by  a 
reference  signal  r(t): 


Wm(s) 


Ym(s) 

R(s) 


which  is  equivalent  to 


ym(t)  =  -  amym(t) +  bmr(t),  or  (5) 

TO)  =—r—R{s)i 

s  +  am 

for  any  bounded  piecewise  continuous  signal  r(t),  where 
am  >  0,  bm  ^  0  are  known.  We  assume  that  am,  bm,  and 
r(f)  are  chosen  so  that  ym(t)  represents  the  desired  output 
signal. 


Before  introducing  the  MRC  law,  we  start  with  two  defini¬ 
tions  and  a  theorem  (proven  in  the  Appendix). 

Definition  1 

A  state  xe  is  said  to  be  an  equilibrium  state  of  the  system 
x  =  /(f,x),  x(f0)  =  x0,  where  x  €  R",/  :  ST  x  B(r)  — > 
Rn,£?  =  [t0, 00),  B(r)  =  {x  S  1"  |  1 1 x 1 1  <  r},  if  f(t,x.e)  = 
0  Vf  >  to-  We  assume  that  f  is  such  that  for  every  xo  €  B(r) 
and  every  t0  £  [0,  00),  the  system  possesses  one  and  only  one 
solution  x(f;to,xo). 

Definition  2 

A  equilibrium  state  xe  is  exponentially  stable  if  there  exits  an 
a  >  0  and  for  every  e  >  0  there  exists  5(e)  >  0,  such  that 


We  will  now  design  a  proper  controller  u(t)  such  that  all 
signals  in  the  closed-loop  system  are  bounded  and  the  system 
output  yp(t)  tracks  the  reference  model  output  ym(t).  The 
control  law  should  be  chosen  so  that  the  closed-loop  plant 
transfer  function  from  the  input  r(t)  to  the  output  yp(t)  is 
equal  to  the  reference  model  transfer  function.  Motivated  by 
this,  we  propose  the  control  law 

C(s ,  0*e)  =  U(s)  =  -klYP(s)  +  k*2R(s)  -  k*3D(s), 

or  equivalently,  in  the  time  domain 


u*(t)  =  - k{yp(t )  +  k2r(t)  -  k3d(t),  (6) 


where  k't,  k2,  k3  are  controller  coefficients  chosen  so  that 


TO 


R(s)  s  +  am 

Eq.  (7)  is  satisfied,  if  we  select 

K  =  ~  (#5  -  am),  k2  = 

which  yields 
1 


TO) 

R(s)  ■ 


(7) 


A/V 


If*  = 

3  P1P2  ’ 


Ps 


( t )  =  TT-TOa  -  am)yp(t)  +  ^TO)  -  yy-^-d(t),  (8) 


P1P2 


P1P2 


P1P2 


provided  of  course  that  Pi,  P2,  P3  0,  i.e.,  the  system 
is  controllable.  Such  a  transfer  function  matching  guarantees 
that  yp(t)  =  ym(t),  Vf  >  t0,  when  yp(t0)  =  ym(to ),  or 
\yP(t)  ~  ym(t)\  ->  0  exponentially  fast  when  ym(to)  ±  2/P(*o) 
for  any  bounded  reference  signal  r(t).  We  also  note  that, 
depending  on  the  parameters  of  some  patients,  this  law  may 
yield  a  negative  control  signal  which  can  not  be  implemented 
in  practice  (corresponds  to  “extraction  of  bivalirudin”  from 
the  patient).  In  such  a  case,  we  need  to  set  a  lower  threshold 
of  zero  for  the  control  signal.  The  final  MRC  control  signal 
becomes  max{0,w(f)}  with  u(t)  defined  as  in  (8). 

We  test  the  performance  of  the  MRC  on  the  data  set  we 
described  in  Section  II-B.  We  only  use  the  test  set  (1/3  of 
the  total)  for  testing  since  the  remaining  training  set  was  used 
for  model  parameter  identification.  As  mentioned  before,  d(t) 
(the  linear  combination  of  physiological  variables  at  time  t ) 
is  a  step-wise  signal  over  time.  By  applying  the  parameter 
identification  method  outlined  in  Sec.  II-B,  we  obtained  both 
population-wide  parameter  values  and  individual  model  pa¬ 
rameter  values. 

We  tested  the  MRC  control  law  on  a  subset  of  patients  and 
the  results  were  qualitatively  the  same  in  each  case.  We  report 
results  from  a  randomly  selected  patient  who  has  identified 
model  parameters  and  available  input  data.  To  that  end,  we 
set  the  reference  parameters  as  am  =  10,  brn  =  700,  r(t)  =  1. 
Choosing  these  values  keeps  the  reference  PTT  value  to  be 
70s,  which  is  in  the  middle  of  the  desirable  range.  We  note  that 
these  parameter  values  are  simply  an  example  and  physicians 
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Fig.  4.  The  effect  of  the  MRC  law  derived  for  and  applied  to  one  randomly  selected  patient. 


have  the  freedom  of  selecting  alternative  values  depending  on 
the  stable  value  and  response  time  they  wish  to  achieve. 

The  effect  of  the  MRC  law  (8)  on  this  randomly  selected 
patient  is  shown  in  Figure  4.  It  can  be  seen  that  driven  by 
inputs  generated  by  the  MRC  law,  the  system  output  quickly 
converges  to  the  reference  output  (top  figure).  The  tracking 
error  (e(t)  =  yp(t)  —  ym{t ))  quickly  converges  to  zero  and 
remains  at  zero  (second  figure).  We  also  obtain  the  control 
signal  which  corresponds  to  the  bivalirudin  infusion  rate  (third 
figure).  The  MRC  control  law  we  introduced  is  robust  to 
the  uncontrollable  signal  d(t)  (bottom  figure).  Although  d(t) 
changes  over  time,  the  control  signal  can  adapt  and  drive  the 
system  to  track  the  reference  signal  closely. 

We  next  evaluate  the  performance  of  the  MRC  scheme  on 
all  patients  in  the  test  set  we  described  earlier.  For  performance 
evaluation,  we  use  three  performance  metrics.  The  first  one  is 
the  Root  Mean  Square  Error  (RMSE),  which  for  patient  i  is 
debited  as 


RMSE,  = 


\ 


1 

Ti 


A* 

t=t\ 


ViniU))2, 


where  t]. ,  tj1  are  the  time  instants  at  which  we  adapt 
the  controller  for  patient  i.  We  debne  RMSE  for  the  whole 
population  of  patients  as  the  average  per  patient  RMSE, 
i.e.,  RMSE=  ^  RMSE,,  where  Nt  is  the  number  of 
patients  in  the  test  set.  We  also  debne  aEi  to  be  the  standard 
deviation  of  the  errors  el(t)  =  ylp(t)  —  ylm{t),  t  =  tj, . . .  ,t?\ 
of  patient  i.  Similarly,  we  debne  07,;  as  the  average  standard 
deviation  of  crEi’s,  i.e.,  aE  =  J2?=i  aEi- 

To  capture  a  notion  of  “relative”  error,  we  also  compute  the 
Normalized  Root  Mean  Square  Error  (NRMSE)  debned  for 
each  patient  i  as 


NRMSE,  = 


N 


1 

Ti 


'52[(ylp(u) 

t=t} 


yin(ti))/ylm(ti)]2- 


As  with  the  RMSE,  we  debne  the  population-wide  NRMSE 
as  the  average  of  NRMSE,; ’s  over  the  patients.  Similarly,  we 
also  debne  cr^Ei  as  the  standard  deviation  of  the  normalized 
tracking  errors  e*(t)  =  ( ylp(t )  -  ylm(t)) / ylm(t)  for  patient  i, 
and  (Tne  as  the  average  of  crNEi’s  over  the  patients. 

Furthermore,  to  illustrate  the  percentage  of  PTT  outliers 
which  are  outside  clinically  safe  bounds,  i.e.,  not  in  the  interval 
[ym(t)  —  10 ,ym{t)  +  10]  during  the  transient  and  not  in  the 
interval  [60s,  80s]  in  steady-state,  the  Risk  Percentage  (RP) 
for  patient  i  is  debned  as 


RPi 


Arrisk 

T 


where  jVJ'*sfc  is  the  number  of  time  instants  t], at 
which  the  PTT  value  of  patient  i  is  outside  the  safe  bounds. 
Then,  RP  is  debned  as  the  average  of  RP,  ’s  over  patients.  We 
note  that  this  metric  is  from  a  clinical  perspective  the  most 
important  in  assessing  the  efficacy  of  our  methods.  Table  I 
reports  the  performance  of  the  MRC  law  for  the  patients  in 
the  test  set. 


TABLE  I 

Performance  of  the  Model  Reference  Control  (test  set) 


value 

RMSE 

0.84 

CTE 

0.82 

NRMSE 

1.20% 

&NE 

1.17% 

RP 

0% 

In  summary,  in  the  case  that  model  parameters  are  known, 
the  MRC  law  tracks  the  reference  signal  quite  well  as  demon¬ 
strated  by  the  low  RMSE  and  NRMSE.  The  RP  value  is  zero, 
which  completely  assures  clinical  safety.  The  corresponding 
standard  deviations  for  RMSE  and  NRMSE  are  small  as  well. 
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Fig.  5.  The  MRC  law  derived  for  one  patient  but  applied  to  another  patient. 


B.  Indirect  Model  Reference  Adaptive  Control  (MRAC) 

As  we  mentioned  in  the  Introduction,  there  is  significant  pa¬ 
tient  variability  in  the  response  to  bivalirudin.  We  have  already 
established  [1],  that  adapting  model  parameters  to  individual 
patients  leads  to  improved  performance.  This  suggests  that  the 
model  structure  is  largely  accurate  but  model  parameters  of  an 
individual  patient  can  deviate  from  population-wide  parameter 
values. 

TABLE  II 

Performance  of  the  Inaccurate  Model  Reference  Control 
(test  SET) 


value 

RMSE 

9.93 

(TE 

3.98 

NRMSE 

14.18% 

&NE 

5.68% 

RP 

61.60% 

To  better  assess  the  effect  of  this  variability,  we  test  the  per¬ 
formance  of  the  MRC  law  derived  using  parameter  values  of  a 
specific  patient  when  applied  to  another  patient  with  different 
model  parameters.  Fig.  5  plots  the  MRC  law  performance  for 
such  a  case.  The  top  figure  shows  that  there  exists  a  large 
gap  between  the  reference  output  signal  and  the  system  output 
signal.  In  addition,  the  system  output  is  outside  the  safe  range. 

We  also  tested  the  MRC  law  derived  using  parameter  values 
of  a  specific  patient  against  all  patients  in  the  test  set.  Table  II 
reports  the  results.  We  note  that  RMSE,  NRMSE,  and  RP  are 


substantially  higher  (and  clinically  unacceptable)  than  those  in 
Table  I.  The  difference  of  course  is  due  to  the  fact  the  results 
of  Table  I  are  obtained  when  using  the  MRC  law  with  the 
correct  model  parameters  for  each  patient,  whereas  Table  II 
results  apply  the  MRC  law  with  parameters  of  some  patient 
to  another  patient.  The  small  values  of  <je  and  <jne  in  Table  II 
indicate  that  performance  of  the  MRC  law  when  used  with  the 
wrong  parameters  is  consistently  poor.  This  situation  should 
obviously  be  avoided  because  overdosing  or  underdosing  is 
very  dangerous  for  the  patients. 

To  address  this  important  issue,  we  next  develop  a  method 
that  first  estimates  the  individual  model  parameters,  and  then 
adopts  the  MRC  law  we  introduced  using  a  certainty  equiva¬ 
lence  principle  [13],  Such  a  control  scheme  is  called  indirect 
Model  Reference  Adaptive  Control  (MRAC)  law. 

By  adding  and  subtracting  —amyp(t)  to  (4),  we  can  obtain 
the  State-Space  Parametric  Model  ( SSPM): 

yP(t )  =  - amyp(t)  +  (arn-/33)yp{t)+l31l32u(t)+f33d(t ).  (9) 

Based  on  (9),  the  series-parallel  estimation  model  [14]  is  given 
by: 

yp(t)  =  - amyp(t )  +  ( am  -  /33 (t))yp(t)  +  pi(tjf}2(t)u(t) 

+  %  {t)d(t),  (10) 

where  yp(t)  is  an  estimated  value  of  yp(t),  and  /3i(f),  P2 (f), 
/33(f)  are  estimates  of  the  system  parameters  [3±,  (32,  and  /33  at 
time  t.  Note  that  in  (10),  yp(t)  is  treated  as  an  input  available 
for  measurement.  By  using  the  certainty  equivalence  principle 
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(cf.  (8)),  we  take  the  control  scheme  structure  to  be: 

u(t)  =  -ki (t)yp(t)  +  k2(t)r(t )  -  k3(t)d(t),  (11) 

where 

h(t)  = 

818)828) 


~  ^ ^ - • 

ihimit) 

In  this  problem,  we  will  estimate  the  product  of  8i(t)  and 
828)  instead  of  estimating  them  separately.  The  model  esti¬ 
mation  error  is  e(t)  =  yp(t)  —  yP(t)  which  implies: 

e(i)  =VP8)  -  Vp{t) 

=  -  ame(t)  +  /33 (t)(yp(t)  -  d(t))  -  812 (12) 

where 

83(f)  =  838)  ~  83,  8128)  =  818)828)  —  8182-  (13) 


k2(t)  = 


83(t) 


818)828) 


We  now  choose  a  Lyapunov-like  function 

V(t)  =  2e(^)2  +  2^^3(*)2  +  2^^12(i)2’ 

which,  with  71,  72  >  0,  is  non-negative  for  all  t.  By  taking 
the  derivative  on  both  sides  of  (14)  we  obtain 


V(t)  =e(t)e(t)  +  —  83(f)  83(f)  +  —8128)8128) 
7i  72 


=  -  am(e(t))2  +  83(f) 


1 


e(t)(yP(t)  -  d(t))  H - /33(f) 


7i 


■  812(f) 


—812(f)  ~  e(t)u(t) 
72 


(15) 


Then,  choosing 


83(f)  =  -li(yP(f)  -  d(t))e(t)  and  812(f)  =  72 e(t)u(t) 

leads  to  V(t)  =  —ame(t)2  <  0.  In  addition,  since  8182  and 
83  are  constants,  (13)  implies  83(f)  =  83(f)  and  812(f)  = 
d(8i(t) 82(f)) / dt.  It  follows  that  we  could  estimate  the  model 
parameters  by: 


838  +  At)  =  83(f)  +  838)^8 
812  (t  +  At)  =  8l2(t) 


d(0  i  («)&(«)) 
dt 


At, 


(16) 


for  small  At.  Then,  we  can  adapt  the  controller  coefficients 
recursively  and  control  the  system  in  real-time  by  using  (cf. 

(ID) 

*/,,  0,m  ~  83(f)  /.\  .  km  .  83(f)  ,/,n  ,,7, 

■u  (t)  = - ^ - yp(t)  +  — - r(t)  —  — - d(t).  (17) 

812(f)  812(f)  812(f) 


Similarly,  as  we  did  earlier  for  the  MRC  law,  we  will  use  the 
control  max{0,  u*(t)}  to  avoid  negative  values. 


Theorem  III.2  Under  the  control  law  (17),  the  tracking  error 
converges  to  0  as  t  — >  00. 

Proof:  By  choosing  such  control  law,  V  (t)  = 

— ame(t )2  <  0,  Vt  >  to-  Since  V(t)  is  bounded  from  below 


and  non-increasing,  it  converges  to  a  constant.  This  implies 
that  —  am  e2(t)dt  =  V"(oo)  —  V(to)  is  bounded,  which  in 

turn  implies  that  e(t)  — >  0  as  t  — >  00  according  to  Barbalat’s 
lemma  [15].  It  also  follows  that  /)3(f),  812(f)  — 1 ►  0  as  t  — >  00. 


One  key  flaw  of  the  adaptive  control  law  (17)  is  that  the 
boundness  of  control  signal  u(t)  can  not  be  established  unless 
we  show  that  k\(t),  k2(t),  fc3(f)  are  all  bounded.  However, 
such  a  control  law  may  generate  estimates  of  8182  arbitrarily 
close  or  even  equal  to  zero,  which  leads  to  the  uncontrollability 
of  the  estimated  model  and  unboundness  of  u(t).  To  avoid 
this  issue,  we  propose  a  modification  to  the  control  law  (17). 
One  method  is  to  modify  the  adaptive  law  for  812(f)  so  that 
adaptation  takes  place  in  a  subset  of  E  which  does  not  include 
the  zero  element.  We  need  to  use  the  a  priori  knowledge  of 
81  8  8li  >  0  and  82  >  8l2b  >  0  to  do  the  projection: 

83(f)  =  -  71  (yP(t)  -  d(t))e(t), 

7 2e(t)u(t),  if  \8i2(t)\  >  8182^ 
or  \8i2(t)\  =  81i82> 
and  e(t)u(t)sga(8i2(f))  >  0, 

0,  otherwise. 


After  modifying  the  adaptive  control  law,  the  time  derivative 
of  the  Lyapunov  function  becomes: 

'  —am(e(t))2 ,  if  \8i2(t)\  >  8182, 
or  \8i2(t)\  =  8[b8ib 

■  =  and  e(t)w(t)sgn(/312)  >  0, 

[  ’  \  -am(e(t))2 

+8i2(t)e(t)u(t),  if  |^i2 (t) |  =  8182 

and  e(f)u(f)sgn(/3i2(f))  <  0. 


Therefore,  V(t)  <  — ame2(t )  <  0,  V<  >  t-o- 

Using  a  similar  argument  as  before,  it  can  be  shown  that  by 
using  this  modified  parameter  estimation  law  (18),  the  tracking 
error  converges  to  zero  driven  by  a  bounded  control  signal. 
Additionally,  we  have  shown  (cf.  Thm.  III.  1 )  that  the  reference 
output  response  is  exponentially  stable,  and  it  follows  that  the 
system  output  can  be  driven  to  the  stable  state  exponentially 
fast. 

We  next  test  the  indirect  MRAC  law  using  the  patient  data. 
The  parameter  values  of  the  reference  model  are  the  same  as 
in  Sec.  III-A.  We  choose  the  population-wide  parameter  values 
83  =  7.9  x  10-4,  and  8182  =  4-22  as  initial  values  of  83(f) 
and  812(f) •  respectively.  The  MRAC  adapts  based  on  these 
estimates  in  real-time.  We  also  set  71  =  72  =  5  x  10-4.  The 
trajectory  of  the  system  under  the  indirect  MRAC  is  shown  in 
Fig.  6. 

Fig.  6  indicates  that  the  system  output  quickly  converges  to 
the  reference  output  and  it  remains  within  the  desired  range 
(top  figure).  The  tracking  error  oscillates  around  zero  (middle 
figure),  but  it  is  not  as  smooth  as  in  Fig.  4.  This  is  due  to 
the  fact  that  the  indirect  MRAC  takes  some  time  to  estimate 
the  system  parameters  first  and  then  adapts  the  controller 
coefficients.  Similarly,  we  can  also  obtain  the  bivalirudin  in¬ 
fusion  rate  (bottom  figure).  Notice  that  although  d(t)  changes 
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llme(days) 


Fig.  6.  The  performance  of  the  indirect  MRAC  law. 


over  time,  the  control  signal  can  drive  the  system  to  track 
the  reference  output  signal  well.  The  performance  of  indirect 
MRAC  for  all  patients  in  the  test  set  is  reported  in  Table  III. 


We  will  devise  a  control  law  that  estimates  these  coefficients 
directly. 

Consider  the  error  derivative: 


TABLE  III 

Performance  of  the  Indirect  Model  Reference  Adaptive 
Control  (test  set) 


value 

RMSE 

5.52 

&E 

2.00 

NRMSE 

7.8S% 

&NE 

2.86% 

RP 

2.30% 

e(t)  =yP{t)  -  ym{t) 

=  -  ame{t)  +  PiP2[-ki(t)yP(t) 

+  k2(t)r(t )  -  k3(t)d(t)], 

where  kft)  =  kft )  —  k*,i  =  1,2,3.  It  follows  that  kft)  = 
ki(t),  i  =  1,2, 3. 

To  design  the  controller,  consider  the  Lyapunov-like  func¬ 
tion: 


We  note  that  in  this  case,  the  values  of  our  three  perfor¬ 
mance  metrics  are  all  higher  than  those  in  the  Table  I.  This 
is  again  explained  by  the  fact  that  adaptation  of  the  indirect 
MRAC  law  is  not  very  fast  due  to  the  parameter  estimation 
step.  Yet,  RMSE,  NRMSE  and  RP  values  are  now  significantly 
less  than  those  in  Table  II  and  they  could  be  considered 
acceptable  in  a  clinical  setting. 

C.  Direct  Model  Reference  Adaptive  Control  MRAC 

In  this  section,  we  focus  on  designing  the  direct  MRAC 
without  estimating  the  model  parameters  first.  Because  the 
individual  model  parameters  are  unknown,  we  can  not  apply 
the  MRC  law  directly.  Based  on  the  structure  of  MRC  law, 
we  propose  an  adaptive  control  law  with  similar  structure: 

u(t)  =  —ki(t)yp{t)  +  k2(t)r[t )  -k3(t)d(t), 

where  ki(t),  k2(t),  and  k3{t)  are  the  estimates  of  MRC 
controller  coefficients  k3,  k'f  and  If  at  time  t,  respectively. 


V(t)  = 


e{tf 


PiP2 


ki{t)2  k2(f)2  k3  (f)2 


271 


272 


273 


By  taking  the  derivative,  we  obtain: 
V(t) 


=  e(t)e(t)  +  P1P2 


ki(t)ki  (t)  k2(t)k2(t)  k3(t)k3(t) 


7i 


72 


73 


=  -ame2(f)  +  f3\f32 
k2(t) 


kl^\ki(t)  -  e(f)7it/p(f)) 


7i 


h(t). 


+  ~^^{k2(t)  +  e(f)72r(f))  H - ~{k3(t)  -  e{t)j 3d{t)) 

72  73 

where  we  have  the  a  priori  knowledge  that  (3 i/32  >  0. 
Choosing 


h(t)  =7i  e{t)yp(t), 
k2{t)  =  -  7 2e{t)r(t), 
k3{t)  =73  e(t)d(t), 
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time(days) 


Fig.  7.  The  performance  of  direct  MRAC  law. 


leads  to  V(t)  =  —ame(t)2  <  0.  Furthermore,  the  controller 
coefficients  can  be  adapted  by 

ki(t  +  At)  =  ki(t)  +  ki(t)At,  i  =  1,2, 3, 

and  our  direct  MRAC  controller  becomes 

u(t)  =  -ki(t)yp(t)  +  k2(t)r(t )  -  k3{i)d(t).  (19) 

As  with  the  previous  two  controllers,  we  will  use  max{0,  u(t)} 
to  avoid  negative  controls. 

We  establish  the  following  result;  we  omit  the  proof  because 
it  is  similar  to  the  proof  of  Theorem  III.2. 

Theorem  III.3  Under  the  control  law  (19),  the  tracking  error 
converges  to  0  as  t  — >  oo. 

We  know  that  the  reference  model  output  is  exponentially 
stable  and  the  tracking  error  converges  to  zero  as  t  increases, 
which  establishes  that  the  system  output  can  be  driven  to 
track  the  reference  output  well.  To  avoid  overdosing  and 
underdosing  risks  that  may  occur  in  the  time  it  takes  the 
MRAC  to  converge  to  the  appropriate  controller  parameters, 
we  can  use  as  initial  estimates  of  these  coefficients  the 
MRC^ coefficients,  namely  Aq(0)  =  2.37,  k2{0)  =  165.75, 
and  k3  (0)  =  1.87  x  10-4,  obtained  from  a  system  model 
parametrized  with  population-wide  parameters. 

We  applied  such  a  direct  MRAC  law  to  the  same  patient 
used  for  generating  Fig.  6.  The  result  is  shown  in  Fig.  7 
(using  71  =  72  =  73  =  0.001).  We  note  that  driven  by 
the  direct  MRAC  law,  the  system  output  quickly  converges 
to  the  reference  output.  Although  the  system  output  oscillates 
around  the  reference  signal,  it  remains  within  the  desired  range 
(top  figure).  The  tracking  error  also  converges  to  zero  (middle 


figure).  We  also  plot  the  bivalirudin  infusion  rate  (bottom 
figure).  Compared  to  the  indirect  MRAC,  the  direct  MRAC 
has  similar  performance  on  controlling  the  PTT.  However, 
the  direct  MRAC  avoids  parameter  estimation  and  estimates 
controller  parameters  directly. 

Table  IV  reports  the  overall  performance  of  the  direct 
MRAC  for  the  patients  in  the  test  set.  Compared  to  the  results 
from  the  indirect  MRAC,  the  direct  MRAC  achieves  lower 
values  for  all  performance  metrics,  i.e.,  RMSE,  NRMSE,  erg, 
one  and  RP.  Notice  the  rather  significant  decrease  of  the  RP 
value,  which,  as  we  argued,  is  clinically  a  top  priority.  This 
table  validates  the  fact  that  the  direct  MRAC  is  a  more  efficient 
and  safer  control  scheme  than  the  indirect  MRAC. 

TABLE  IV 

Performance  of  the  Direct  Model  Reference  Adaptive 
Control  (test  set) 


value 

RMSE 

0.80 

(7E 

0.80 

NRMSE 

1.15% 

one 

1.14% 

RP 

0.09% 

IV.  Conclusions 

Based  on  a  specific  dynamic  system  model  of  bivalirudin 
acting  in  cardiac  surgical  patients,  we  developed  two  meth¬ 
ods  for  synthesizing  a  controller  to  regulate  the  bivalirudin 
infusion  rate  and  induce  a  PTT  within  a  desirable  range.  The 
first  method  assumes  that  the  model  parameters  are  available 
and  develops  a  control  law  that  tracks  a  physician  specified 
reference  output  signal.  Our  second  method  considers  patients 
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for  which  past  clinical  records  are  sparse  and  accurate  model 
parameters  are  not  readily  available.  It  develops  an  indirect 
control  scheme  (indirect  MRAC)  that  first  estimates  the  model 
parameters  and  then  adapts  the  corresponding  controller  based 
on  these  estimates.  Alternatively,  a  direct  control  scheme 
(direct  MRAC)  that  adapts  the  controller  without  estimating 
the  model  parameters  first  is  also  developed.  Testing  of  these 
schemes  against  actual  patient  data  from  a  hospital,  shows  that 
the  direct  MRAC  is  more  efficient  than  the  indirect  version. 

The  methods  we  developed  can  be  seen  as  key  steps  towards 
automation  of  dosage  decisions  in  a  hospital  setting,  which 
can  help  eliminate  errors  and  neutralize  the  inexperience  of 
residents  who  are  currently  responsible  for  these  decisions. 
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Appendix 


The  solution  to  Eq.  (5)  is 


t 

ym{t)  =  $(t,t0)ym(t0)  +  j  T>(f,r)6mr(r)dT, 
to 


where  <£>((,  t)  =  e~a,n<-t~r'>  is  the  state  transition  function 
in  this  problem.  Since  r(t )  =  Cr,  which  is  a  constant,  the 
solution  to  (5)  can  be  written  as 


ym(t o)  - 


h  C 

um.K-^r 


b  rn  Cr 


(20) 


Equation  (20)  indicates  that  ffu  Z/m^o))  — >  bmCr  which 
is  a  constant,  as  t  — ►  oo.  In  addition,  using  Definition  1,  it  can 
be  easily  verified  that  yme  =  b m Gr  is  the  equilibrium  state  of 
our  reference  system.  Furthermore,  \ym(t;  f0, ym(to))~  Vme\  = 
|g-am(t— ; to)  {ym(t0)-yme)\  =  |t/m(f0)-2/me|e-am(*-to),  Vt  > 
to-  Therefore,  by  Definition  2,  it  follows  that  the  reference 
model  equilibrium  state  yme  is  exponentially  stable. 
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